Application of a multi-site mean-field theory 
to the disordered Bose-Hubbard model 
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We present a multi-site formulation of mean-field theory applied to the disordered Bose-Hubbard 
model. In this approach the lattice is partitioned into clusters, each isolated cluster being treated 
exactly, with inter-cluster hopping being treated approximately. The theory allows for the possibility 
of a different superfluid order parameter at every site in the lattice, such as what has been used in 
previously published site-decoupled mean-field theories, but a multi-site formulation also allows for 
the inclusion of spatial correlations allowing us, e.g., to calculate the correlation length (over the 
length scale of each cluster). We present our numerical results for a two-dimensional system. This 
theory is shown to produce a phase diagram in which the stability of the Mott insulator phase is 
larger than that predicted by site-decoupled single-site mean-field theory. Two different methods 
are given for the identification of the Bose glass-to-superfluid transition, one an approximation 
based on the behaviour of the condensate fraction, and one of which relies on obtaining the spatial 
variation of the order parameter correlation. The relation of our results to a recent proposal that 
both transitions are non self-averaging is discussed. 



I. INTRODUCTION 



The experimental realization of a gas of ultracold atoms in an optical lattice [l| has provided a unique opportunity 
to exploit a remarkably close connection between experimental and theoretical studies of interacting Bose condensed 
systems. In one sense, one may think of the experiments as constituting a "quantum simulator" of strongly interacting 
particles on a lattice Q described by the Hubbard model Q. With the subsequent addition of speckle-generated 
disorder Q to these experiments, or by the generation of multichromatic lattices Q, a window into the world of 
strongly interacting disordered physics was opened. 

The seminal paper of Fisher, et al. @, in 1989 predates the production of Bose-condensed "cold" atomic gases by 
several years. The focus of that paper is the disordered Bose-Hubbard model, including the thermodynamically stable 
quantum phases, the properties of these phases, and the transitions between them. Many theory papers have followed 
since this work, the number of such papers is growing (nearly) exponentially since the realization of this model via 
cold atom experiments. 

One of the simplest and most often employed methods to analyze a model such as the disordered Bose-Hubbard 
^ model is mean- field theory. The paper of Fisher, et al. [6|, used such a theory to provide the first qualitative phase 
diagram of this system. Many mean-field theory investigations, some of which we now mention, have followed this 
■ paper. Our contribution, presented in this report, is the application of a so-called multi-site mean-field theory to 
investigate some of the same questions addressed in Ref. [6|, with the ultimate goal of making comparisons between 
. such a theory and experiment. As we outline below, the hoped-for improvement of such a multi-site theory is the 
i—\ ' incorporation of some spatial correlations that are neglected in conventional single-site site-decoupled mean-field 
analyses. 

As stated above, the present work is an example of a multi-site mean-field theory, and to put this approach into 
perspective it is important to first discuss the set of results found from the familiar site-decoupled single-site mean-field 
theory. The original paper of Fisher, et al. [(|, provided the qualitative aspects of the phase boundaries separating 
the relevant quantum phases, namely the Mott insulating, Bose glass, and superfluid phases. Some time ago, using an 
RPA-based approach, this phase diagram was obtained quantitatively Q • Following the many experimental advances 
discussed above, this theory was re-examined by several groups. We mention first the spatially-averaged single-site 
mean- field theory of Krutitsky, et al. Q. As discussed in the next section, this formulation allows one to obtain 
the phase boundaries analytically. However, since one is treating each site as being equivalent to all other sites, the 
resulting theory does not include any of the effects of coherent back-scattering that give rise to localized states, and 
correlations between the order parameter at neighbouring sites are completely ignored (see Eq. below). This 
limitation may be overcome by allowing for a different order parameter at every site, something that was achieved 
in a series of detailed papers (below we refer specifically to only three papers in this series, Refs. [9TJIlj|. since these 
are most relevant for comparisons to our results). These latter authors mapped one aspect of the problem onto an 
effective non- interacting Anderson localization problem and thereby found a simple analytical device to obtain 
(numerically) the Mott insulator-Bose glass phase very accurately, albeit within single-site theory. Then, to obtain 
the Bose glass-to-superfluid transition they analyzed the superfluid stiffness. They distinguished the Bose glass phase 



2 



from that of the superfluid by stating that the condensate fraction, f c , is non-zero in both of these phases, while the 
superfiuid stiffness is non-zero only in the superfluid region. However, unlike the enormous lattices that were used by 
these authors to obtain the Mott insulator-to-Bose glass transition lines, when examining the Bose glass and superfluid 
phases one must solve a set of self consistency equations for each site on a lattice, and this limits the sizes of the 
systems that can be studied. Therefore, one is faced with the difficulties of completing calculations for many lattice 
sizes, thereby allowing for a finite-size scaling analysis to be completed. Also, one must average over many different 
complexions of disorder. Unlike the statements in their paper Q, when we completed comparative single-site numerics 
in the Bose glass and superfluid phases, we found considerable finite size effects, as well as substantial variations with 
differing complexions of disorder. As we discuss below, a proposal exists that could explain this behaviour. 

The potential pitfalls of finite-size scaling and averages over complexions of disorder are avoided in the very elegant 
and creative theory of Bissbort and Hofstatter, so-called stochastic mean-field theory [l2|, [l3|. These authors have 
found that for (and, so far, only for) single-site mean-field theory one may perform exactly (albeit numerically) the 
avera ging over different complexions of disorder in the thermodynamic limit. Unfortunately, as discussed in detail in 
Refs. |l3l fl3|. the stochastic mean-field theory method is one that can only examine the thermodynamic limit. These 
authors argue, in direct contrast to the work of Refs. [H, [ll|, that the condensate fraction in the Bose glass phase 
is identically zero, a result intimately tied to the existence of rare Lifshitz regions. Part of the motivation for our 
work is a comparison to the finite systems studied experimentally, for example the observed Bose glass phase in three 
dimensions [15]], as it is important to be able to understand large but finite systems. 

We have completed large-scale numerical studies of both the single-site approach of Refs. [H, [ll| as well as the 
stochastic mean-field theory of Refs. [12, EH, and have found differences between these theories. Therefore, this 
provided further motivation to go beyond single-site theory to see if the correct physics could be identified. 

There is also a comparison of the single-site mean-field theory with Quantum Monte Carlo studies available in the 
literature [l6l Il7j . Comparing both position and momentum distributions of particles (including the presence of a 
harmonic trapping potential), it has been shown that for a three dimensional Bose-Hubbard model the mean-field 
approximation results closely resemble the exact solution when the superfluid fraction is non-zero. It is suggested 
that the lack of particle-particle correlations in the single-site formulation of the mean-field theory is a deficiency that 
seems to affect primarily Mott insulator and Bose glass phases. 

The idea of extending single-site theories to multiple sites, or clusters, is not, of course, original. In a sense the 
same idea has been used, e.g., when the very successful DMFT (dynamical mean-field theory) [181 ] was broadened to 
cellular or extended DMFT [l9|,|20j]. That is, such work represents an attempt to include more spatial correlations 
into a theory that ignores the wave-vector dependence of the self-energy. The same cluster extension has also been 
adapted to DMFT for disordered systems - see the paper of Potthoff and Balzar [2l| in which the Luttinger self-energy 
functional expansion technique is employed. Also, a similar idea was applied to the diatomic Bose-Hubbard model, 
as we outline below. Therefore, the inclusion of such effects is a natural extrapolation from site-decoupled mean-field 
theories that have been applied in other circumstances before. 

So, as discussed, the familiar limitations of mean-field theory apply to site-decoupled mean-field theory, notably 
the absence of spatial correlations. For example, for the Bose-Hubbard model one can't calculate pair correlations of 
the order parameter at differing sites as this corresponds to simply the one-body density matrix. Our contribution is 
the extension of site-decoupled mean-field theory for the disordered Bose-Hubbard model to a cluster, or equivalently 
multi-site mean-field theory. Our theory treats exactly the spatial correlations on each cluster (each cluster has open 
boundary conditions), thereby giving us access to spatial correlations over length scales up to the largest distances 
spanning a cluster. In one dimension this allows for us to have clusters up to a size of six (eight for homogeneous 
lattices), whereas for two dimensions even though one can examine clusters up to a size of 3x3 for homogeneous 
clusters [22j], for spatially disordered systems we were required to work with clusters of size 2x2, although one can 
also work with clusters of 3x2 sites using the technique of basis truncation [26j . which is described later in this paper. 
Like the site-decoupled single-site mean-field theory discussed above, in our multi-site work for disordered systems we 
also allow for a different (self-consistent) order parameter at every site in the lattice. 

Our paper is organized as follows. In the next section we present the formalism associated with multi-site mean-field 
theory applied to the disordered Bose-Hubbard model. In this section we also include the formal details associated 
with the spatially-averaged variant of multi-site mean-field theory, a natural extension of Ref. Q . In the following 
section we present our numerical results of applying this theory to the two-dimensional disordered Bose-Hubbard 
model, including an extrapolation of the spatially-averaged variant of single-site mean-field theory of Ref. [8] to a 
multi-site formulation. This is followed in the last section by a discussion of our results, their relation to previously 
published data, and further extensions that are possible with this theory. 
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II. FORMALISM 

The Hamiltonian under study is the disordered Bose-Hubbard model, which is defined (in the grand canonical 
ensemble) as 

t = n- fiti ' = -t^2 + H - c ) + ( £i ~ /i ) hi + \ H ~ W 

where i,j label the sites of the optical lattice, &i and at annihilate and create a boson on site i, respectively, 
denotes that i and j are nearest neighbours, and the number operator for site i is denoted by hi — a|aj. The 
near-neighbour hopping frequency is denoted by t, and all hopping beyond nearest neighbours is suppressed. The on- 
site (Hubbard) interaction energy is denoted by U, and consistent with recent experiments (23[, interactions between 
bosons on different sites are considered to be small and are ignored. The disorder is introduced into the model 
considered here via the on-site energies E{, and we focus on a uniform distribution of energies from the range of — W/2 
to +W/2. The chemical potential fi controls boson number. In what follows below, we restrict our considerations to 
zero temperature. 



A. Site-decoupled mean-field theory: 

Our work examines the application of a multi-site mean-field theory to this model. To understand the differ- 
ences with the familiar and oft-employed site-decoupling mean- field theory (e.g., see Ref. Q), note that this latter 
approximation can be derived by writing 

hi = fa + S&i where fa = (*o|oi|*o) (2) 

with j^o) representing the ground state, and the (spatially homogeneous) superfluid order parameter being denoted 
by fa Then, ignoring terms in the hopping Hamiltonian that are second order in 5a, the resulting model Hamiltonian 
(denoting site-decoupling by sd) for a homogeneous lattice, where all = and therefore the order parameter is 
independent of lattice site, namely fa — fa is given by 

jC*d = J2tCf (3) 

i 

where 

JCf = —hi(hi-l)- iihi - zt (hi + af) (4) 

where z is the coordination number for the lattice under consideration. Equation ((3|) makes clear that this formulation 
leads to a sum over single-site Hamiltonians. Further, this means that the form of the ground-state wave function is 

i*o>=ni^> ( 5 ) 

i 

When one considers a disordered system, the order parameters will, in general, all be different (except in the Mott 
insulating regions, for which the order parameter fa is zero everywhere). Then one must generalize Eq. Q to be 

= ^fiiih, - 1) + [si - fj)hi -t(a,i + af) ^ 4>i+s (6) 

where the last sum is over 5 a sum over all sites within one near- neighbour distance to site i. Therefore, one again 
obtains a Hamiltonian that is site decoupled [1C s i d depends on creation and annihilation operators only for site i), and 
a ground-state wave function of the form of Eq. ([5]). However, such a formulation does include spatial correlations in 
the sense that for a given complexion of disorder (that is, for a given set of on-site energies £j), one must solve the 
self consistency equation at every site. That is, using Eq. ([5]) the site-dependent order parameter is found from 

fa = (\& |fit|*o) = (ipi\ai\ipi) (7) 

and as is seen from Eq. (6), the ground state of each ICf d depends on the order parameters on neighbouring sites. 
Last, we note the well-known result that the minimum value of the grand potential (at T = 0) 

n sd = (* |/C sd |* ) (8) 
is achieved when the self-consistency conditions given in Eq. ([7]) are satisfied. 
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B. Multi-site mean-field formulation: 



The multi-site version of mean-field theory for the homogeneous Bose-Hubbard model has been stated previously 
in the literature [25j | . and these authors used this formulation to examine a one-dimensional diatomic lattice. They 
obtained data by analyzing one diatomic "cluster" and examined the stability of so-called loop-hole phases. In a 
subsequent paper these authors stated that the convergence with cluster size (see below) was poor, and therefore 
preferred to examine their problems using different approaches. The work of Mcintosh, Zaremba and one of the 
present authors [22[ showed how one obtains improved phase boundaries when one uses a multi-site formulation, in 
particular when examining two and three-dimensional lattices. Further, in the analysis of the one-dimensional diatomic 
lattice the stability conditions on various "loop-hole" , or "fractional lobe" , phases and mass-density wave states were 
seen to be changed when more than one diatomic cluster was considered. Below we discuss the extrapolation of these 
ideas to disordered systems. 

The simplest derivation of this formulation can be given in one dimension, and for clarity we detail this formalism for 
such situations. Consider a one-dimensional chain with N s lattice sites, the sites being labelled by i = 0, 1, . . . , N s — 1, 
with periodic boundary conditions. Let the chain be decomposed into N c clusters with L sites per cluster, and 
N a = N C L. Label the clusters by I = 0, 1, . . . , N c — 1, and note that any given lattice site can now be indexed by 
i = £ + IL with £ = 0, 1, 1. Then, write the (exact) grand Hamiltonian as 



(9) 



where JCj is given by 



^■i ~ J2e=o....,L-i 



U 



nt+i L {nt+iL - 1) + {ee+iL - m) n i+ 



IL 



-t Y,t=o,...,L-2 {al +IL ae+i+iL + a\ +1+IL a l+IL 



(10) 



Each K,°j is the exact grand Hamiltonian for a cluster of length L with open boundary conditions. Therefore, one should 
think of Eq. ([9]) as representing a sum of single-cluster Hamiltonians (JC°j) that are then coupled via hopping terms 
between the near-neighbour clusters (Vj,j). One can think of the inter-cluster interactions as being perturbations on 
the sum over cluster Hamiltonians. 

The idea behind the multi-site mean-field approximation is to apply mean-field decouplings to the inter-cluster 
hopping terms, V/,j. That is, applying Eq. @ only to these terms, and denoting the order parameters at each end 
of the cluster by </>i, where i = L — 1 + IL or i = L + IL for the left and right ends of the cluster, one finds that the 
multi-site mean-field approximation (denoted by the superscript ms) for the inter-cluster interaction term is given by 



V? 



ms 
1+1 



-t 



a-L+iL 



l L+IL 



l+IL 



O-L-l+IL 



it 

l L-l+IL 



?L+IL 



Then, the multi-site mean-field Hamiltonian is seen to be given by 



(ii) 



(12) 



The main difference between the site-decoupled and multi-site mean-field theories will be made clear below. For 
now, simply note the differing forms of the ground-state wave functions. For the site-decoupled formulation one 
obtains the form shown in Eq. (JS|), whereas for the multi-site formulation one has 



iM> )=niv>/> 



(13) 



where \tpx) is an energy eigenstate of the cluster Hamiltonian K,j. That is, the eigenstates of ICj include the spatial 
correlations present in a fully interacting cluster of length L with open boundary conditions, and specifically one 
sees that these eigenstates depend on the hopping parameter, t. If instead one treats the hopping as a perturbation 
in the site-decoupled formulation, the unperturbed eigenstates do not depend on the hopping frequency. This begs 
the question, does the inclusion of some effects of hopping into the zeroth-order eigenvectors lead to more accurate 
estimates of the phase boundaries and quantum phases? This was, in part, our motivation for this study. 

The geometry associated with the generalization of this idea to higher dimensions, and in particular to the focus of 
this paper, namely the two-dimensional square lattice, is straightforward - the site and cluster labels are now vectors, 
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and an example partitioning of the lattice is shown in Fig. [T] The only technical complication is that every site on 
the "edge" of a cluster has its own order parameter, and for a disordered system all such order parameters are, in 
general, different. 

An important difference is seen immediately - the zcroth-ordcr cluster Hamiltonian for a site-decoupled theory does 
not depend on the hopping frequency, t, and does not include any spatial correlations between the (disordered) on-site 
energies. 




FIG. 1: A portion of a two-dimensional square lattice that is partitioned into 2x2 clusters - only the central cluster is shown 
in its entirety. The dashed lines show the division of the lattice into such clusters. In each cluster there are four sites, labelled 
by A,B,C,D, and each cluster is denoted by a subscript, which runs from 1 to 9. 



The form of the cluster states that are eigenvectors of fC°j can be written as 

n max n max Ti max 

w = EE- E ^o,n 1 ,...,n i _ 1 (4o) no (4 1 ) ni ---(4,_ 1 r i - i io) 



(14) 



no— Oni— riL-i—0 



and numerically one must choose some finite n max that leads to fully converged results. We have found that for the 
region of the phase diagram associated with the first Mott lobe, namely < p, < 1, n max = 3 is sufficient when 
t < 0.08 - all of our results shown in this paper are within this range. 

The dimensionality of the Hilbert space grows very rapidly with increasing cluster size and so does the computational 
time. Therefore, in order to go beyond clusters of 2x2 sites, one needs to reduce (truncate) the number of basis states 
even further. The challenge is to decide which states are important and should be retained and which can be discarded. 
Of course, one still wants the results to be converged. Following the work of Schmitt, et al. we have adopted a 
truncation scheme which allows for a physically motivated and controlled approach that discards those states which 
contribute the least to the particular cluster's ground state. As a simple measure for the importance of individual 
number states one can use the expectation value of the Hamiltonian. The key idea is to discard only those basis states 
for which the expectation value of the Hamiltonian is greater than some value. Those states have a considerably 
smaller contribution to the final ground state of a particular cluster, and this procedure leaves us with a given 
percentage of the initial basis set. We found that for clusters of 3x2 sites it is enough to retain only 30% of the "full" 
basis set, "full" corresponding to n max = 3 of Eq. (fl"4)) (as discussed in the previous paragraph). Having in mind that 
the truncation scheme works well [26| in the regions where the ground state is composed of few dominant basis states, 
one has to remember that the results for the superfluid phase will be less reliable. (One needs the full real space basis 
set to accurately describe the superfluid phase, since it is a phenomenon observed in the momentum space.) 



C. Quantum phases: 



We will work in the strong-coupling limit in that we assume that U is the largest energy scale, and instead of 
working with K, and its mean-field variants we work with the dimensionless 

= - (is) 
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and also express the various material parameters as scaled energies, namely 

~ t 11 Si 

* -> * = — fi -> A = jj £i -> £i = ^ (16) 

Letting the (scaled) energy VF characterize the strength of the disorder, the phase diagram for this model can be 
obtained in (t, jl, W) space. To begin, we present a brief discussion of the phases that appear in such phase diagrams. 

For a homogeneous lattice with only near-neighbour hopping and on-site repulsion, one finds that Mott insulating 
phases are found to correspond to an integral number of bosons per site. The superfluid phase is then found in the 
(i, ji) plane to correspond to (f> 7^ 0. For non-homogeneous lattices the situation is more complicated. Specifically, 
when = one finds that the Mott insulator phase is an incompressible phase that corresponds to integral densities 
of bosons per lattice site. For ordered non-homogeneous lattices, for example for a diatomic lattice in one dimension, 
one finds a variety of density- wave phases [22|, |27J, (so that the average density per site is integral, but the local 
bosonic densities need not be integral) that are also incompressible and therefore Mott insulators. In addition to 
such behaviour, even for binary disorder one finds Mott insulating phases but now of a different nature [28j ; these 
incompressible phases do not even have an average boson density per site that is integral. Our analysis focuses on 
non-homogeneous lattices with a disorder that is not binary (from now on in our discussion it is implicit that we 
are not addressing the situation of binary disorder) , and the transition to non-Mott insulating phases is qualitatively 
different. 

Formally, one notes that the exact grand Hamiltonian of Eq. (fTJ) commutes with the total number operator, and 
therefore in the absence of degeneracies or symmetry-breaking terms being added to Eq. (fTJ), the eigenstates are 
also number eigenstates. For systems with spatially disordered on-site energies, £j, this is also true of both the site- 
decoupled (Eq. ([3])) and multi-site (Eq. <\V2^i ) mean- field Hamiltonians when the order parameters fa are zero on all 
sites. However, when the mean-field terms are included they act as a source of symmetry-breaking interactions, and 
now the ground state is found to not correspond to number eigenstates (cither for a single site or for a multi-site 
ground states - see Eqs. (|5I13[) ). Therefore, the resulting phases are compressible. Also, as seen from Eq. (8), in such 
a situation one obtains self-consistent solutions for which fa ^ 0. 

For a homogeneous lattice the compressible phase is known to be a superfluid. However, for spatially disordered 
systems it is known that between the Mott insulating and superfluid phases a new quantum state arises, a phase 
called the Bose glass phase. Several important results concerning this transition were summarized in the introduction 
to this paper. Here we present the methodology that we have employed to determine the phase boundary separating 
the Mott insulating and Bose glass phases in (i, jl, W) space. 

To identify the phase boundary of the Mott insulating and superfluid phases for the homogeneous system one may 
use a variety of methods. For example, one may use perturbation theory 29] on a Landau-type expansion of the 
grand potential fl in the order parameter fa Simply, when the coefficient of the fa 2 term is found to change from 
being positive to negative (for say constant jl with increasing t), the order parameter becoming non-zero signifies a 
continuous phase transition from the Mott insulator to superfluid phases. This is the method that Mcintosh, Zaremba 
and one of us have used [22| in analyzing various ordered systems. 

A variant of this approach has been used for disordered systems by treating the site-decoupled mean-field theory 
with a spatial average over disorder Q. To be concrete, if one approximates the order parameter to be spatially 
homogeneous for the disordered system, namely fa — > fa, and ignores the correlations between neighbouring sites, so 

that fa(f>j — > (0) , then using the same approach as for a homogeneous system (for a given distribution of on-site 
energies) in a site-decoupled mean-field approximation one may derive [§[ (in fact, analytically) the loss of stability 
of the Mott insulating phase with increasing disorder. Below we explain how this approach can be extended to 
our multi-site mean-field formulation; such results serve as the most simplistic approach in using our formalism to 
determine the locations of the phase boundaries. 

The approach that utilizes the full spatial structure of the on-site energies and order parameters that we have used 
can be explained as follows. Again for clarity, consider a one-dimensional lattice with L sites per cluster. For a given 
complexion of disorder any lattice will have different values of the order parameter at every site (outside of the = 
Mott insulator region). Then the phase boundary in such an approach can be derived by consideration of a stability 
matrix. That is, beginning with 

i {I, J) 

= ^^? + f6 
i 
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where (for simplicity with the notation, restricting our attention to one-dimensional chains) 



&L+IL 



l L+IL 



vl-i+il 



O-L-X+IL 



h-l+IL j <PL+IL 



Then, representing the decoupled (that is, 6 = 0) eigenstates via 



i 



(18) 



(19) 



with all eigenstates of Eq. (fT9l) being labelled by the index a, and Eq, |^q) being the ground-state energy and 

eigenvector of the decoupled system respectively, one determines the ground state of JC ms to lowest order in the 
perturbation 6, namely 



i*o> = i*s> 



\*P) 



where 



(20) 



(21) 



The self-consistency equation, to lowest order in 6, can then be written as 
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(22) 



since one notes, from Eq. (|18p . that varies linearly with 0,. This quantity, C, is what we refer to as the stability 
matrix. That is, when the magnitude of any eigenvalue of C exceeds one, the Mott-insulating phase, corresponding 
to 0, = Vi, is found to be unstable. 

For the homogeneous lattice one finds that the Mott insulating phase is unstable to the formation of a superfiuid. 
However, for the disordered Bose-Hubbard model one finds an intervening Bose glass phase. As is well known [(I [30|. 
the appearance of the condensate in an interacting system may be identified by a non-zero condensate fraction, denoted 
by f c , which may be evaluated based on the largest eigenvalue of the one-body density matrix in the thermodynamic 
limit. This latter quantity is defined as follows. First, let 2V& denote the total number of bosons, namely 



N b = J2 n i 



(* |ni|* ) 



(23) 



Then, the one-body density matrix may be defined by 



p\ 3 



(24) 



Now consider the multi-site mean-field theory for a one-dimensional crystal, and, to be concrete, consider a cluster 



size of L = 2. Noting that when the order parameters 4 
same form as in Eq. (|13[) . it is then simple to show that 



1 



are non-zero the ground state wave function is still of the 
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Note that for the homogeneous system all order parameters are equal, and one has that Ni, = nN s 
largest eigenvalue of this matrix would then be given by 



1 



(«o»i) 
nN s 



(Ns ~ 2)0 2 
nN s 



Therefore the 



(26) 
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and thus 

f c = lim A max = — (27) 
jv s — >oo n 

as expected. (Identical results are obtained for any homogeneous system treated with multi-site mean-field theory, 
with any cluster size.) For the disordered system, the largest eigenvalue can be extracted numerically, and finite-size 
scaling should allow for one to obtain f c . The scaling form follows from Eq. (|26[) . and below and in our figures we 
denote this finite-size corrected condensate fraction to be denoted by f c . 

The quantum phase that appears at the Mott-insulating phase boundary is known to be the Bose glass phase [f|. 
The condensate fraction in this phase, at least for finite systems, is non-zero. However, unlike the homogeneous 
system, for a disordered system the superfluid stiffness, p s , may be zero while the condensate fraction, f c , is non-zero. 
In the thermodynamic limit this is purported to be associated with islands of condensate that are not phase coherent 
with one another [l3j . Below we show how the complete phase diagram can be (approximately) determined from the 
condensate fraction alone. 



III. RESULTS FOR THE 2D SQUARE LATTICE WITH BOX DISORDER 

We now turn to the evaluation of the phase boundaries and an examination of the properties of the quantum phases 
for a two-dimensional square lattice having spinless bosons described by the disordered Bosc-Hubbard model, given 
in Eq. ([T]), and present our numerical results for both the site-decoupled single-site mean- field theory and for our 
multi-site formulation. We have solved this problem for a variety of lattice sizes such that the number of lattice sites 
times the number of complexions of disorder is 40,000, and the lattice sizes range up to 200x200 sites. We discuss 
the "convergence" that we obtain with these choices below. In our multi-site work we use mainly clusters of size 2x2, 
since the dimensionality of the Hilbert space of a cluster grows very rapidly with increasing cluster size. To reduce the 
number of basis states for homogeneous lattices one can use group theory to find equivalent lattice sites, and hence 
equivalent basis states [22J. However, in the presence of disorder the C4 point group symmetry is broken, so no such 
formal device is available. Therefore, to go beyond clusters of 2x2 sites we have used a basis truncation scheme that 
was introduced elsewhere pril |. and which was described earlier in the previous section of the present work. Through 
this approach we have been able to extend our study to examine the results that follow from the use of clusters of 3 x 2 
sites. Our results are limited to the region of the phase diagram associated with the first Mott lobe, namely < jl < 1, 
when determining the phase boundaries (see Eq. (|22|)^ and self-consistent solutions for the order parameter at each 
site of a given lattice. Also, we present results for the case of a uniform (or box) distribution of on-site energies with 
W — 0.5 (see Eq. (IT51) '). More details on the state labelling (and related numerical issues) can be found in |3l| . 

The Mott insulator-to-Bose glass phase boundary obtained within the chosen range of jl is given in Fig. [51 For 
the single site mean-field theory (see Fig. [2^,) we show three different lattice sizes - 25 complexions of disorder for a 
40x40 lattice, 4 complexions of a 100x100 lattice and a single complexion of a 200x200 lattice. Additionally, for the 
200x200 system, we show the phase boundary of a different realization of disorder, but for only three values of the 
chemical potential; the chemical potentials used are those for which we saw the 100x100 lattice to have the biggest 
variations in the phase boundary results. It is clear from this figure that not only does one obtain a range of phase 
boundaries at smaller lattice sizes, but also for much larger lattices. In contrast to the claim of Ref. who studied 
many lattices and complexions of disorder in one dimension and found no dependence on either of these (and claimed 
that a chain of length 100 did not suffer from finite-size effects), and two dimensions in other published work [l(J, in 
our own numerical study of site-decoupled single-site mean-field theory we find a strong dependence of these phase 
boundaries on lattice size and complexions of disorder. (However, for the condensate fraction for large systems our 
single-site theory numerics are consistent with their statements - see below.) 

We show the results for our multi-site formulation in Fig. [2}d, with the same number of complexions of disorder 
for each of the lattice sizes studied for single-site theory. For all lattice sizes and all complexions of disorder we find 
that the area of the phase diagram corresponding to the Mott insulator phase is greater when calculated using our 
multi-site theory than that predicted by single-site theory, a result consistent with studies on homogeneous lattices 
0. 

If one examines Fig. [5] at some fixed chemical potential but with an increasing hopping frequency, in particular 
for chemical potentials in the lower half of the Mott lobe (namely around 0.3 < ji < 0.45), one finds a considerable 
variation of the location of the phase transition even for our largest lattices. An identical spread of the phase boundaries 
is found for single-site theory, even for lattices with 200x200 sites. This result, that one does not find a convergence 
to a single phase transition [321 ] . is not consistent the self averaging that one might expect for large lattices. Therefore, 
for most experimental systems (read: finite with less than 40,000 lattice sites), for a given complexion of disorder one 
can expect some variation in the stability of the Mott insulating and Bose glass phases, depending on both the extent 
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40x40 □ 200x200 

• 100x100 ■■■■ Spatially Averaged 

FIG. 2: Mott insulator/Bose glass phase boundaries of a two-dimensional (a) site-decoupled single-site mean-field theory 
(SSMFT) and (b) multi-site mean-field theory (MSMFT), where the clusters are 2x2 sites, for a box distribution of disorder 
of strength W = 0.5. The results shown are for lattice sizes of 40x40, 100x100 and 200x200, with 25, 4 and 1 complexion 
of disorder, respectively, and are represented using the line types/point styles given in the legend. For example, we are 
superimposing 25 stability curves for the 40x40 lattices. The spatially-averaged phase boundary is found using the approach 
of @] for single-site mean-field theory, and as described in the main body of the paper for 2x2 clusters. 

of the harmonic trap (thus dictating the (approximate) effective size of the lattice), and on the particular complexion 
of disorder that is produced in a given experiment. An interesting aspect of this observation is that it is seen to be 
consistent with a recent renormalization group analysis of this transition [33| - for the disorder fixed points, they 
predict the absence of self averaging in this model. 

In addition to the numerical results of Fig. [2J one may extend the spatially averaged disorder theory of Krutitsky, 
et al. 8], to the multi-site formulation, as we now explain. Referring to the above-outlined perturbation theory, if 
one wishes to determine the variation of the energy (more specifically, the grand potential) when the effects of the 
perturbation are included (see Eq. p7|) ). as in Refs. [1, [H|, one obtains a Landau-type thermodynamic potential for 
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each cluster of the form 

flj = + Aij<fti(f)j + higher-order terms (28) 

The expansion coefficients depend on the various on-site energies, e, of Eq. (fT6|) . and paralleling the approach of 
Ref. Q we approximate that for all i,j 

J^~^ 2 (29) 

where the spatially averaged superfluid order parameter is used to replace all local order parameters appearing in 
Eq. (|28l) . Therefore, the Mott insulator-to-Bose glass transition is found to correspond to the crossing of zero (say 
with increasing i at some fixed jl) of the second-order coefficient. (Of course, Eq. (|2"9"]l is satisfied exactly in the Mott 
insulator.) For a homogeneous lattice this formalism is presented elsewhere [22] . and for a disordered lattice using 2x2 
clusters the extrapolation of Q to our multi-site formulation requires that one average over all possible on-site energies 
on each of the four sites, something that can only be completed numerically. The single-site result of @ is plotted in 
Fig. [5^,, and for the multi-site theory the resulting phase boundary is shown as a dotted line in Fig. ^p. Notably, it is 
seen to be much closer to the numerical results described above for large lattices than is obtained with site-decoupled 
single-site mean-field theory. That is, by including the spatial correlations within each multi-site cluster one obtains 
a reasonably accurate representation of the MSMFT phase boundary. (Unlike the single-site theory Q , we are not 
able to complete the required integrals analytically and therefore can not write down a simple analytical expression 
for the phase boundaries of a spatially averaged multi-site theory.) Clearly, with increased computing power this idea 
can be extended to larger clusters, thereby providing one with a much better disorder-averaged phase boundary, if 

A O 

the Kj eigenstates can be calculated. 

The variation of the condensate fraction with hopping frequency, correcting for finite-size effects, as described at 
the end of the previous section, is shown in Fig. [3] for two values of the chemical potential; one below the first Mott 
lobe (which for W = 0.5 is found to be bounded by 0.25 < ji < 0.75), namely jl = 0.1, and one inside, that is 
pL = 0.4. Again, both single- and multi-site mean-field theory results are shown for a variety of lattice sizes for which 
the number of lattice sites times the number of complexions of disorder is 40,000. We find that in both single- and 
multi-site formulations of the mean-field theory the condensate fraction seems to converge to a fixed value independent 
of the complexion of disorder when large enough lattices are used, but only when this fraction is away from zero; this 
behaviour is similar to what was found in the numerics of Ref. [9|. Also, if one examines these curves closely in the 
vicinity of the Mott insulator-to-Bose glass phase transition, when the condensate fraction is very close to zero (for 
jl = 0.4), one finds that the condensate fraction becomes nonzero precisely at the values of the hopping frequency 
found from the instability theory (see Fig. [5]). That is, when f c first becomes non-zero for a given lattice size varies 
with complexion of disorder, and does not converge to a fixed value. Below we refer to this hopping value as t c \. 
(For completeness, note that i c \ = for ji = 0.1.) As is clear from these plots, we find that the implementation of 
a multi-site theory results in qualitatively similar changes to those of the phase boundaries that are found from the 
single-site formulation. That is, the onset of f c (say f c greater than ~ 0.01) is shifted to higher values of t, similar to 
the shift of MI/BG phase boundary discussed above - see Fig. [2j 

The above results have led us to identifying a simplistic approximation that shows how one can estimate not only 
when the Mott insulator loses stability (as discussed above, this corresponds to t c i), but how one can also obtain 
an estimate of the phase boundary separating the Bose glass and superfluid phases, considering only the variation of 
the condensate fraction. That is, at larger hopping frequencies the disorder is found to not suppress phase coherence 
throughout the lattice (see below), and one obtains a superfluid. We shall identify the latter transition by t c2 . 

To make clear how our approach works, in Fig. 0] we show the variation of this quantity for values of the chemical 
potential for which the Mott insulator phase does not exist (t c \ — 0), viz. for jl = 0.1, and also for which it is stable 
only below some critical value of t c \ > 0, viz. for jl — 0.4. One can estimate the critical values t®| an d t™, at which 
the system transforms from being in the Bose glass phase to that of a superfluid, as follows. Recall from above that 
one may expand the grand potential in terms of a spatially homogeneous order parameter, (j>, phenomenologically as 
a Landau-type expansion, namely 

n = n° + l^(t c -i)* ® 2 + |(0) 4 (30) 

From this it follows that when the order parameter first becomes non-zero, one has that 



(31) 
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- 20x20 — 100x100 

40x40 ■■■■ SMFT 



FIG. 3: The variation of the condensate fraction, f c , with hopping frequency, t, for W = 0.5 and (a) ji — 0.1 and (b) jl — 0.4 - 
the appropriate finite-size corrections for the condensate fraction are included. For single-site mean-field theory (SSMFT) the 
results for lattices of, respectively, 20x20, 40x40, and 100x100 sites with 100, 25 and 4 complexions of disorder are shown. The 
legend shows linetypes/colours to which the particular lattice size results correspond. One sees that (consistent with Ref. [9|]) 
there is virtually no change in f c for sufficiently large lattices - the 40x40 curves are nearly coincident, and the 100x100 curve 
lies on top of them - providing that the condensate fraction is not close to being zero. Note that our results from stochastic 
mean-field theory (SMFT), denoted by the dotted line, agree with the large lattice results for the single-site theory, but only 
for the condensate fraction f c > 0.1 (see inset). In addition, we show the results obtained from our multi-site mean-field theory 
(MSMFT) for clusters of 2x2 sites. The lattice sizes and numbers of complexions of disorder are the same as studied for 
single-site theory. 

a common mean-field result. Then noting that for a homogeneous system the superfluid fraction, f s , is equal to the 
condensate fraction, our previous result of f c ~ <f> 2 shows that one then obtains that the superfluid fraction first 
develops linearly in the hopping frequency. From Fig. 2] one sees that when the condensate fraction first becomes 
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FIG. 4: The variation of the condensate fraction, f c , for both SSMFT and MSMFT, with hopping frequency, t, for W = 0.5 
and (a) p, = 0.1 and (b) fi = 0.4 for one complexion of disorder. Both the SSMFT and (2x2) MSMFT results have been 
obtained on a 100x100 lattice for the same complexion of disorder. The results for (3x2) MSMFT have been obtained on a 
30x30 sites lattice with 30% basis truncation (see text for more details). The linear fit of the data provides an estimate of the 
Bose glass-to-superfluid phase transition, an d i respectively, for SSMFT and MSMFT. The Mott insulator-to-Bose glass 
transition at tjf and t™^ are also shown for fi — 0.4; for jl = 0.1 these i c i quantities are both zero. 



sizable it grows linearly with t, and therefore one can find a simplistic approximation for t C 2 by extrapolating back 
from the linear region. The results of this procedure are shown in the figure and are denoted with (red) dotted lines. 
The results shown refer to single-site mean-field theory, t*2, and multi-site mean- field theory with clusters of 2x2 
sites, denoted by t^Jf . 

To lend support to the validity of this approach, and to demonstrate another advantage of using MSMFT, we have 
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FIG. 5: The variation of the correlation length £ (from Eq. (|33[0 with hopping frequency t, for one particular complexion of 
box disorder of strength W = 0.5, evaluated using our multi-site mean-field theory with clusters of 3x2 sites. The results have 
been obtained on a 30x30 sites lattice with 30% basis truncation (see text for more details). The error bars show asymptotic 
standard error of the correlation length fit parameter. Note that the peak of the correlation length for jl = 0.4 agrees with the 
location of t C 2 w 0.043 estimated from the linear extrapolation of the previous figure (see the linear extrapolation of the 3x2 
cluster data in the previous figure). 



examined and quantified the spatial correlations present in our multi-site ground state wave functions - of course, 
this information is not available when a site-decoupled single-site theory is employed. That is, we have calculated 

C(\T i -T j \) = {S$\S$ j ) 1 = (32) 

where the overline corresponds to an average over all pairs of sites. Then this function is fit to a functional form 
appropriate when the system is in the vicinity of a phase transition, namely 

ccfc-^.aHs^za (33 , 

(Note that for our 2x2 clusters one can access only the distances of 1 and y/2 lattice constants. However, one has 
many bonds at each of these distances.) The case of 2x2 clusters is particularly easy for the analysis, since then the 
correlation length is given by £ = — l/lnC(l). We find that for jl — 0.4 the correlation length takes its maximum 
value near the hopping frequency corresponding to i^ s , found from Fig. 21 consistent with our simplistic approach of 
identifying the Bose glass-to-superfhiid phase transition shown in that figure. 

We have also evaluated f c and the corresponding correlation function using larger clusters, namely of size 3x2 sites, 
and our results are shown in Fig. [SJ The implementation of this larger lattice gave us two additional distances of 2 and 
\fh lattice constants. Because the computational time grows very rapidly with increasing cluster size, this analysis has 
been limited only to lattices of 30x30 sites. We found that the parameters describing asymptotical behaviour of the 
correlation function (see Eq. (|33p) closely resemble the ones found for 2x2 clusters. Just as for the smaller clusters, 
the maximum of the correlation length is found to be consistent with the linear fit of the f c data for jl = 0.4, as 
discussed in the caption. For jl = 0.1 the comparison of the phase transition identification between these two methods 
is not as good, and is possibly due to the fact the variation of the correlation length with hopping is very broad and 
not peaked. Therefore, while we caution that our simplistic approximation for identifying the Bose glass-to-superfluid 
phase transition by a linear extrapolation of the condensate fraction data should be only treated as an approximation, 
our results on the correlation length, at least for p, = 0.4, lend support for its validity. Furthermore, if larger clusters 
can be handled numerically, then the use of a MSMFT and the evaluation of the pair correlation function could 
provide new information about the behaviour of the phase coherence in the Bose glass phase, a topic that we discuss 
more below. 

As mentioned earlier, our numerical result, and previous work [H, [llj], on finite lattices, that the condensate fraction 
is non-zero in the Bose glass phase, is in contrast to the stochastic mean-field theory, which applies to systems in the 
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thermodynamic limit, as espoused in Refs. [H, [Hj]. Our results in Figs. 13141 notably the inset of Fig. [3] for p, = 0.1, 
make clear that the curvature of f c vs. t in the immediate neighbourhood of does not change as one increases 
system size. In stochastic mean-field theory this rounding is changed to a discontinuity of the slope with f c reduced 
to zero below the critical value t c i- As mentioned above, we believe that one can address this apparent disagreement 
through an analysis of the spatial correlations of the order parameter [3~ij , something that requires a multi-site theory 
beyond what is contained in Fig. [SJ but such a resolution of this difference would require numerical results for much 
larger clusters. 

In Fig. [5] we show histograms and position dependencies of both the order parameters (a-d) and site's occupation 
numbers (e-h) in the Bose glass (left panel) and superfluid (right panel) phases obtained from the single- and multi-site 
formulations of mean-field theory. Single-site results are shown in figures a, b, e, f, and multi-site results in figures 
c, d, g, h. For completeness, we also show the random on-site energies Ei for the analyzed 100x100 lattice (Fig. O). 
Shades of grey reflect the magnitude of the analyzed quantities (for values, see gradation scales next to each figure). 
We choose the values of the hopping integral to correspond to condensate fractions of f c « 10 -4 for the Bose glass 
phase, and to f c w 0.05 for the superfluid phase. 

We find that the maps of the order parameters (Fig. [5b, d), as well as the site's occupation numbers (Fig. Hf,h), in 
the superfluid phase for the single-site and multi-site mean-field theories are quite similar. This suggests the benefits 
of using the simpler single-site theory when studying the superfluid phase (apart from the inability of this formulation 
to give information about the correlation length). 

The situation is more interesting in the Bose glass phase. In Figs. |B^,,c we show the position dependencies of the 
order parameters. From these maps one can see that even though for the majority of sites the order parameters are 
equal to zero, some sites can still develop a non vanishing order parameter. (Of course, it is this behaviour which 
leads to the non-zero condensate fraction in the Bose glass phase.) The sites for which this quantity is smaller than 
the convergence criterion (in our numerics, this corresponds to sites with <\>i < 10~ 4 ) have been denoted as blank 
(white). The results are similar for both single- and multi-site formulations, although the multi-site formulation leads 
to order parameters which change less abruptly (smoother edges), something that is not unexpected. In contrast to 
this, the local boson densities obtained within the single-site mean-field theory are almost uniformly distributed (and 
equal to one) over all lattice sites, whereas for the multi-site theory the densities show much greater variation. Again, 
this is what one would expect. For example, when we look at the results from the multi-site mean-field theory we see 
that the inclusion of intra-cluster hopping leads to nonuniform occupation of all sites (in general) for isolated clusters, 
since, in the areas where the order parameter has not been developed, bosons can only hop between sites belonging 
to a single cluster while inter-cluster hopping is suppressed. To minimize the ground state energy particles need to 
reorganize themselves, and since the on-site energy, Si, differs from site to site, so does the particle density. In the 
areas where the order parameter is non-zero bosons can hop between clusters, and this allows them to move outside 
of regions with high £j (e.g., examine the results corresponding to site (32,56) in Fig. [BJ). For the single-site theory 
bosons can move around the lattice only when the order parameter is non-zero, and this leads to uniform distribution 
of particles where (f>i = 0. In the superfluid phase the inter-cluster hopping is allowed on the entire lattice (the order 
parameters are always non-zero). Therefore the inhomogeneity of the local bosons densities is even stronger - to 
minimize the ground state energy, bosons are no longer restrained to particular cluster and can hop around the entire 
lattice. Nevertheless, there are areas where the order parameter is small, and these areas for the single-site mean-field 
theory show a much weaker variation of the particle density than for the multi-site formulation, simply because the 
effective hopping probability, t<pi, is suppressed. This is reflected in a much narrower local boson density histogram 
for SSMFT than for MSMFT. 

For completeness, in the same figure we also show the order parameter probability distribution obtained within the 
stochastic (single-site) mean-field theory, discussed in the introduction [T3 - tl4| - this data corresponds to the dotted 
lines in Fig. |6Jd, normalized to allow for a direct comparison of these results with SSMFT. We found that even though 
the distributions agree deep inside the superfluid region of the phase diagram (for example, for fl = 0.4 one needs the 
hopping frequency of t > 0.045), in the vicinity of the Bose glass-to-superfluid phase transition they are qualitatively 
different. The standard deviation of probability distribution obtained within the stochastic mean- field theory is much 
smaller, in this region, than the one we get from the single-site mean-field theory solved on a finite lattice. This is 
expected, since within the stochastic theory the order parameter vanishes in the Bose glass phase. 

In Fig. [7] we show the momentum distribution of bosons, in all three different phases, given by 

^ = ^E elk ' (r, ~ rj) <^> ( 34 ) 

The results shown have been obtained for a single complexion of disorder on a 100x100 lattice. (This is the same 
realization of disorder potential as the one used for the position space results in Fig. [5]) We found that in the Mott 
insulator phase (Fig.[7k,) the momentum distribution obtained within the single-site theory is uniform - all states with 
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FIG. 6: Histograms and position dependencies of the boson densities and order parameters obtained on a lattice of 100x100 
sites. In panels (a) and (e) results from single-site mean-field theory (SSMFT) are shown, and in panels (c) and (g) results 
from multi-site mean-field theory (MSMFT) are shown, all for the system existing in the Bose glass phase. Similarly, in panels 
(b) and (f) SSMFT results, and in panels (d) and (h) results for MSMFT, are shown for the system being in the superfluid 
phase. The distribution of the on-site energies for this particular complexion of disorder is presented in the lower figure (i). 
The results have been evaluated for W = 0.5 and p, = 0.4, where the hopping frequencies are chosen such that in the Bose glass 
phase f c « 10 -4 , and in the superfluid phase f c w 0.05. The hopping frequencies, corresponding to those condensate fractions, 
take the following values: (a,e) £=0.03721; (b,f) t= 0.04004; (c,g) i = 0.04128; (d,h) i = 0.04397. Note that in the Bose glass 
phase on the majority of the lattice sites the order parameter is vanishingly small, smaller than the convergence criterion, and 
those regions have been denoted with blank (white) colour. The normalized results for the <f> probability density function from 
stochastic mean- field theory (dotted line) are also shown in (a,b). 
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FIG. 7: Momentum distribution of bosons in the first Brillouin zone, along X-F-M path, obtained from the single- and multi- 
site mean-field theories for W = 0.5 and jl — 0.4 on a 100x100 lattice. Figure (a) shows the results in the Mott insulator 
phase, where t — 0.035; (b) in the Bose glass phase where t — 0.03721 for the single-site and t = 0.04128 for the multi-site 
formulation; (c) in the superfluid phase for t — 0.04004 and t — 0.04397, respectively. (These choices are discussed in the text.) 
The solid/dashed lines denote the average over corresponding directions of the momentum wave- vector and the dots show 
actual data points for SSMFT/MSMFT. In (c) the original data points, before averaging, are shown, where (blue/red) open 
circles/squares refer to SSMFT/MSMFT. Also, in (c) it should be noted that the k = data is truncated, since it corresponds 
to no ~ 475 — 480, and as discussed in the text this corresponds to the macroscopic occupation of the ground state in the 
superfluid phase. (Note the reduction of for the SSMFT from 1 to 0.95 away from k = 0, and a similar reduction for 
MSMFT.) 



momenta k are equally occupied. Note that in this region of the phase diagram all order parameters are identically 
equal to zero, which leads to (a\aj) = {hi)5ij. Therefore, we find that nt = Nb/N s — n, which for the set of t, jl and W 
parameters that we chose, ends up to be n = 1. However, we observe that the momentum distribution resulting from 
the multi-site formulation has a characteristic bell shape that is present throughout the entire phase diagram, even in 
the Bose glass and superfluid phases. This is a consequence of including the intra-cluster hopping in the evaluation 
of the eigenstates of the zeroth order Hamiltonian (see Eq. (fTU)) ), which, as emphasized above, leads to the multi-site 
theory retaining some of the inter-particle spatial correlation effects. Because of this, the momentum distribution of 
bosons in the Mott insulator phase is not a simple sum of the occupation numbers on each site, instead containing 
other terms for which the sum in Eq. (|34[) runs over the individual cluster sites. These effects are completely washed 
out if the single-site formulation is used. In the Bose glass phase (Fig. [7b) the disorder starts to play an important 
role, as bosons can scatter off of the disorder potential, which results in varying particle numbers being present around 
the k = momentum. As one expects, the presence of disorder breaks all symmetries of this distribution; however, 
the inversion symmetry stays intact (not shown). A further increase of the hopping frequency causes the system to go 
into the superfluid phase, and therefore the ground state starts to be macroscopically occupied, which is reflected in 
the momentum distribution by a peak near the k = momentum (Fig.[7J:). In the figure, this leads to no being around 
475-480 (truncated in the figure for clarity), and which is compensated by the reduction of rik at other momenta such 
that the boson densities sums are similar in all panels. 

Taken together, the above results makes clear that it is in the Bose glass phase that the advantages of using MSMFT 
instead of SSMFT are most noticeable. However, as seen in the last figure, and in our discussion of the correlation 
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length, there are quantities for which the MSMFT provides new information. 

IV. DISCUSSION 

We have derived a multi-site mean-field theory for the disordered Bose-Hubbard model, and compared the results 
found from such an approach with those obtained from the more familiar site-decoupled single-site theory. The main 
differences between the single-site and multi-site formulations of mean-field theory are: (I) One includes the effects of 
intra-cluster hopping, thereby producing a perturbation theory that includes intra-cluster particle-particle correlation 
effects in the unperturbed eigenstates (see Eq. (JTOj) ) . (II) One incorporates some aspects of the spatial variation of the 
on-site energies in the unperturbed eigenstates (the cluster eigenstates of Eq. (|19p depend on intra-cluster complexions 
of disorder), while in single-site theory for the unperturbed eigenstates there is only one on-site energy. 

We found that the multi-site theory predicts a stable Mott insulator phase for greater hopping frequencies than 
the single-site formulation. This is consistent with parallel studies of the similar problem for homogeneous lattices 
[22} . Additionally, for both SSMFT and MSMFT theories we showed that the stability of the Mott insulator phase 
depends on the complexion of disorder - even for lattices of size 200 x 200 sites our results did not converge to a single 
value of t c i, the Mott insulator to Bose glass transition. It would be interesting to see if this dependence on the 
complexion of disorder is also observed in the experimental systems. (Note that it is straightforward to include the 
effects of a harmonic trapping potential into our MSMFT numerics.) Of course, we are limited to finite lattices, so no 
rare Lifshitz regions can be included, and other theoretical approaches are more appropriately suited for examining 
their role in determining the phase diagram of the system in the thermodynamic limit. This limit may also lead to 
the results of stochastic MFT theory being correct for the condensate fraction, namely that it vanishes in the Bose 
glass phase, although such a conclusion would need to be shown to be robust under the extrapolation of this inventive 
approach to a multi-site mean-field theory. 

We found that in the Bose glass phase the position distribution of particles obtained within multi-site theory is 
correlated with the underlying random potential of the on-site energies, whereas the single-site theory predicts that 
the vast majority of sites are uniformly occupied with a given number of bosons per site (only sites which developed 
a non-zero order parameter have varying occupation numbers). This is clearly a result of incorporating the spatial 
variation of the on-site energies in the eigenstates of each isolated cluster. 

We have also shown that these two formulations of the mean-field theory give rise to two qualitatively different 
momentum distributions of particles. The multi-site formulation predicts a characteristic bell shape thereof, a feature 
that is clearly a product of the inclusion of the intra-cluster hopping in the zero order perturbation Hamiltonian (see 
Eq. (|34p ). The single-site theory, however, does not show such a behaviour since it does not incorporate any of the 
particle-particle correlations effects. 

We also presented a way of identifying - approximately - the Bose glass-to-superfluid phase transition and therefore 
completing full phase diagram of the Bose-Hubbard model, solved on a finite size lattice for a given complexion of 
disorder, from the condensate fraction alone. By examining the correlation length we have seen that this approximation 
is expected to be accurate near the tip of Mott lobe, but less accurate away from it. 

Taken together, these results suggest that the MSMFT better accounts for the variations of the boson density and 
order parameter distribution that are expected in the Bose glass phase, and it seems here that the potential utility of 
this approach could best be exploited. Future work will discuss the dynamics that are found, in the Bose glass phase, 
using MSMFT. The differences between the two formulations in the superfluid phase do not seem to be as dramatic, 
although quantities such as the correlation length and the momentum distribution are accessible via MSMFT. 
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